import numpy as np
from math import pi
import scipy
from scipy import signal
from matplotlib import pyplot as plt
from scipy.io.wavfile import write
f1 = 900 #formant frequency
b1 = 200 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
r = np.exp(-b1*pi*T)
r
theta = 2*pi*f1*T
theta
plt.rcParams['figure.dpi'] = 600
plt.figure(figsize=(12, 12))
plt.plot(r*np.cos(theta), r*np.sin(theta), 'bx')
plt.plot(r*np.cos(-theta), r*np.sin(-theta), 'bx')
angle = np.linspace( 0 , 2 * np.pi , 150 )
radius = 1
x = radius * np.cos( angle )
y = radius * np.sin( angle )
plt.plot(x, y, 'r', alpha = 0.2)
plt.axvline(0, c = 'r', ls='--', alpha=0.2)
plt.axhline(0, c = 'r', ls='--', alpha=0.2)
plt.plot(0, 0, 'ro')
plt.legend(["The poles of H(z)"])
def getFilterDenCoef(F1, B1, Fs):
r, theta = np.exp(-b1*pi*T), 2*pi*f1*T
return np.array([1, -2*r*np.cos(theta), r*r])
num = np.array([1])
den = getFilterDenCoef(f1, b1, fs)
s = "{a0} + {a1} z^-1 + {a2} z^-2".format(a0 = den[0], a1 = den[1], a2 = den[2])
print("H(z) = \n")
print(" "*(len(s)//2) + "1\n" + "-"*len(s)+"\n"+s)
w, h = scipy.signal.freqz(b=num, a=den, fs = fs)
plt.figure(figsize=(15, 8))
plt.plot(w, 20*np.log(np.sqrt(h*np.conjugate(h))))
plt.xlabel("Frequency")
plt.ylabel("dB(Magnitude)")
plt.title("Magnitude plot of H(z) dB vs Frequency")
plt.grid(True)
w, mag, phase = scipy.signal.dbode((num, den, 1))
plt.figure()
plt.plot(w*fs/(2*pi), mag) # Bode magnitude plot
plt.xlabel("Frequency")
plt.ylabel("Amplitude in dB")
plt.title("Magnitude response of the system")
plt.figure()
plt.plot(w*fs/(2*pi), phase) # Bode phase plot
plt.xlabel("Frequency")
plt.ylabel("Phase(Degrees)")
plt.title("Phase response of the system")
plt.show()
def computeImpRes(ini, samples, r, theta):
#computing the impulse response from the difference equations
y = np.zeros((samples, 1))
y[0], y[1], y[2] = ini
for i in range(3, samples):
y[i] = 2*r*np.cos(theta)*y[i-1] - r*r*y[i-2]
return y
yImp = computeImpRes([1, 2*r*np.cos(theta), 4*r*r*np.cos(theta)*np.cos(theta) - r*r], 200, r, theta)
plt.figure(figsize=(15, 12))
t = [i/fs for i in range(len(yImp))]
plt.plot(t, yImp, '.')
plt.axhline(0, c = 'r', alpha=0.3)
plt.xlabel("Time(seconds)")
plt.ylabel("Amplitude of signal")
plt.title("Impulse response in time domain")
f1 = 900 #formant frequency
b1 = 200 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
r = np.exp(-b1*pi*T)
theta = 2*pi*f1*T
F0 = 160
dur = 0.02
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0 ##approximating location of the peaks of the triangular impulse train
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
def computeResponse(x, samples, r, theta):
#computing the impulse response from the difference equations
assert x.shape[0]==samples
y = np.zeros((samples, 1))
y[0] = x[0] #assumed casual
y[1] = x[1] + 2*r*np.cos(theta)*y[0]
y[2] = x[2] + 2*r*np.cos(theta)*y[1] - r*r*y[0]
for i in range(3, samples):
y[i] = x[i] + 2*r*np.cos(theta)*y[i-1] - r*r*y[i-2]
return y
def triagInput(P, samples, width = 4, start = 0):
#start = {0, 1} whether to start at 0 or the next non-zero sample
x = np.zeros((samples, 1))
peaks= []
for i in range(start, samples+start):
if i%P==0:
x[i] = 1 #populating the peaks
peaks.append(i)
for j in peaks:
for k in range(width):
if j-k>=0:
x[j-k] = 1 - (k/width)
if j+k<samples:
x[j+k] = 1 - (k/width)
return x
plt.subplot(1, 2, 1)
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, computeResponse(triagInput(P, samp, width = width), samp, r, theta))
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.plot(timeSpace, computeResponse(triagInput(P, samp, width = 4), samp, r, theta))
plt.axvline(1/F0, linestyle = '--', alpha = 0.3, linewidth = 1)
plt.legend(["Triangular impulse train", "System output of source-filter"],fontsize=6)
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Output of source-filter system")
F0 = 160
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y = computeResponse(triagInput(P, samp, width = 4), samp, r, theta)#computing the output
scaled = np.int16(y/np.max(np.abs(y)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn2.wav', fs, scaled)
Please listen to qn2.wav file for the audio
Waveforms : We see F0 manifest as the pitch of the sound, that is, the entire 1/F0 = 0.00625 seconds duration window repeats F0 times in a second. (We have a triangular pulse now so components of the next pulse manifest in the previous hence we see that the output signal rises above a littl ebefore this!) F1 is the frequency of the oscillations in this window. We see that is is decaying and this decay is due to our impulse repsonse.
We hear a constant low pitch 'a' sound
f1 = 300 #formant frequency
b1 = 100 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 120
dur = 0.02
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r = np.exp(-b1*pi*T)
theta = 2*pi*f1*T
plt.subplot(1, 2, 1)
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, computeResponse(triagInput(P, samp, width = width), samp, r, theta))
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.plot(timeSpace, computeResponse(triagInput(P, samp, width = 4), samp, r, theta))
plt.axvline(1/F0, linestyle = '--', alpha = 0.3, linewidth = 1)
plt.legend(["Triangular impulse train", "System output of source-filter"],fontsize=6)
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Output of source-filter system")
F0 = 120
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y = computeResponse(triagInput(P, samp, width = 4), samp, r, theta)#computing the output
scaled = np.int16(y/np.max(np.abs(y)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn3a.wav', fs, scaled)
f1 = 1100 #formant frequency
b1 = 200 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 120
dur = 0.02
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r = np.exp(-b1*pi*T)
theta = 2*pi*f1*T
plt.subplot(1, 2, 1)
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, computeResponse(triagInput(P, samp, width = width), samp, r, theta))
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.plot(timeSpace, computeResponse(triagInput(P, samp, width = 4), samp, r, theta))
plt.legend(["Triangular impulse train", "System output of source-filter"],fontsize=6)
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Output of source-filter system")
F0 = 120
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y = computeResponse(triagInput(P, samp, width = 4), samp, r, theta)#computing the output
scaled = np.int16(y/np.max(np.abs(y)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn3b.wav', fs, scaled)
f1 = 300 #formant frequency
b1 = 100 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 180
dur = 0.02
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r = np.exp(-b1*pi*T)
theta = 2*pi*f1*T
plt.subplot(1, 2, 1)
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, computeResponse(triagInput(P, samp, width = width), samp, r, theta))
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.plot(timeSpace, computeResponse(triagInput(P, samp, width = 4), samp, r, theta))
plt.legend(["Triangular impulse train", "System output of source-filter"],fontsize=6)
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Output of source-filter system")
F0 = 180
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y = computeResponse(triagInput(P, samp, width = 4), samp, r, theta)#computing the output
scaled = np.int16(y/np.max(np.abs(y)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn3c.wav', fs, scaled)
Please listen to qn3a.wav for (a) part, qn3b.wav file for the (b) part and qn3c.wav file for the (c) part
f1 = 730 #formant frequency
f2 = 1090
f3 = 2440
b1 = 100 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 120
dur = 0.01
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r1 = np.exp(-b1*pi*T)
theta1 = 2*pi*f1*T
r2 = np.exp(-b1*pi*T)
theta2 = 2*pi*f2*T
r3 = np.exp(-b1*pi*T)
theta3 = 2*pi*f3*T
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)
plt.subplot(1, 2, 1)
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, y3)
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
F0 = 120
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)#computing the output
scaled = np.int16(y3/np.max(np.abs(y3)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn4_a_120Hz.wav', fs, scaled)
f1 = 270 #formant frequency
f2 = 2290
f3 = 3010
b1 = 100 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 120
dur = 0.01
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r1 = np.exp(-b1*pi*T)
theta1 = 2*pi*f1*T
r2 = np.exp(-b1*pi*T)
theta2 = 2*pi*f2*T
r3 = np.exp(-b1*pi*T)
theta3 = 2*pi*f3*T
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)
plt.subplot(1, 2, 1)
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, y3)
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
F0 = 120
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)#computing the output
scaled = np.int16(y3/np.max(np.abs(y3)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn4_i_120Hz.wav', fs, scaled)
f1 = 300 #formant frequency
f2 = 870
f3 = 2240
b1 = 100 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 120
dur = 0.01
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r1 = np.exp(-b1*pi*T)
theta1 = 2*pi*f1*T
r2 = np.exp(-b1*pi*T)
theta2 = 2*pi*f2*T
r3 = np.exp(-b1*pi*T)
theta3 = 2*pi*f3*T
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)
plt.subplot(1, 2, 1)
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, y3)
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
F0 = 120
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)#computing the output
scaled = np.int16(y3/np.max(np.abs(y3)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn4_u_120Hz.wav', fs, scaled)
Now, we go for the high pitch tones
f1 = 730 #formant frequency
f2 = 1090
f3 = 2440
b1 = 100 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 220
dur = 0.01
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r1 = np.exp(-b1*pi*T)
theta1 = 2*pi*f1*T
r2 = np.exp(-b1*pi*T)
theta2 = 2*pi*f2*T
r3 = np.exp(-b1*pi*T)
theta3 = 2*pi*f3*T
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)
plt.subplot(1, 2, 1)
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, y3)
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
F0 = 220
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)#computing the output
scaled = np.int16(y3/np.max(np.abs(y3)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn4_a_220Hz.wav', fs, scaled)
f1 = 270 #formant frequency
f2 = 2290
f3 = 3010
b1 = 100 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 220
dur = 0.02
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r1 = np.exp(-b1*pi*T)
theta1 = 2*pi*f1*T
r2 = np.exp(-b1*pi*T)
theta2 = 2*pi*f2*T
r3 = np.exp(-b1*pi*T)
theta3 = 2*pi*f3*T
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)
plt.subplot(1, 2, 1)
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, y3)
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
F0 = 220
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)#computing the output
scaled = np.int16(y3/np.max(np.abs(y3)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn4_i_220Hz.wav', fs, scaled)
f1 = 300 #formant frequency
f2 = 870
f3 = 2240
b1 = 100 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 220
dur = 0.01
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r1 = np.exp(-b1*pi*T)
theta1 = 2*pi*f1*T
r2 = np.exp(-b1*pi*T)
theta2 = 2*pi*f2*T
r3 = np.exp(-b1*pi*T)
theta3 = 2*pi*f3*T
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)
plt.subplot(1, 2, 1)
plt.plot(timeSpace, triagInput(P, samp, width = width))
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, y3)
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
F0 = 220
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)#computing the output
scaled = np.int16(y3/np.max(np.abs(y3)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn4_u_220Hz.wav', fs, scaled)
Please listen to the following files for the sounds
/a/ @120Hz - qn4_a_120Hz.wav
/i/ @120Hz - qn4_i_120Hz.wav
/u/ @120Hz - qn4_u_120Hz.wav
/a/ @220Hz - qn4_a_220Hz.wav
/i/ @220Hz - qn4_i_220Hz.wav
/u/ @220Hz - qn4_u_220Hz.wav
def glottalPulse(input, a):
y = np.zeros_like(input)
y[0] = input[0]
y[1] = input[1] + 2*a*y[0]
y[2] = input[2] + 2*a*y[1] - (a*a*y[0])
for i in range(3, len(input)):
y[i] = input[i] + 2*a*y[i-1] - (a*a*y[i-2])
return y
We do the above for sound 'a' @120Hz
f1 = 730 #formant frequency
f2 = 1090
f3 = 2440
b1 = 100 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 120
dur = 0.02
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r1 = np.exp(-b1*pi*T)
theta1 = 2*pi*f1*T
r2 = np.exp(-b1*pi*T)
theta2 = 2*pi*f2*T
r3 = np.exp(-b1*pi*T)
theta3 = 2*pi*f3*T
yI = glottalPulse(triagInput(P, samp, width = width), 0.9)
y1 = computeResponse(yI, samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)
plt.subplot(1, 2, 1)
plt.plot(timeSpace, yI)
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train after being shaped by the glottas")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, y3)
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
F0 = 120
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)#computing the output
scaled = np.int16(y3/np.max(np.abs(y3)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn4_a_120Hz_with_glottal.wav', fs, scaled)
Please listen to qn4_a_120Hz_with_glottal.wav
def GlottalLip(input, a,b ):
y = np.zeros_like(input)
y[0] = input[0]
y[1] = input[1] - b*input[0] + 2*a*y[0]
y[2] = input[2] - b*input[1]+ 2*a*y[1] - (a*a*y[0])
for i in range(3, len(input)):
y[i] = input[i] - b*input[i-1] + 2*a*y[i-1] - (a*a*y[i-2])
return y
f1 = 730 #formant frequency
f2 = 1090
f3 = 2440
b1 = 100 #bandwith
fs = 16000 #sampling frequency
T = 1.0/fs #samples
F0 = 120
dur = 0.02
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
r1 = np.exp(-b1*pi*T)
theta1 = 2*pi*f1*T
r2 = np.exp(-b1*pi*T)
theta2 = 2*pi*f2*T
r3 = np.exp(-b1*pi*T)
theta3 = 2*pi*f3*T
yI = GlottalLip(triagInput(P, samp, width = width), 0.9, 0.8)
y1 = computeResponse(yI, samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)
plt.subplot(1, 2, 1)
plt.plot(timeSpace, yI)
plt.xlabel("Time in seconds")
plt.ylabel("Amplitude")
plt.title("Triangular impulse train after glottal-lip filter")
plt.subplot(1, 2, 2)
plt.plot(timeSpace, y3)
plt.xlabel("Time in seconds")
plt.title("System output of source-filter")
F0 = 120
dur = 0.5
timeSpace = np.linspace(0, dur, int(dur*fs))
P = fs//F0
samp = len(timeSpace)
width = 2#no of samples on either side of the peak for the triangular input
y1 = computeResponse(triagInput(P, samp, width = width), samp, r1, theta1)
y2 = computeResponse(y1, samp, r2, theta2)
y3 = computeResponse(y1, samp, r3, theta3)#computing the output
scaled = np.int16(y3/np.max(np.abs(y3)) * 32767)#scaling the output to 16bit integers for writing to wav file
write('qn4_a_120Hz_with_glottallip.wav', fs, scaled)
Please listen to qn4_a_120Hz_with_glottallip.wav
Adding jitter to the sound
freqDelta = np.random.rand(samp)